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ABSTRACT 

We report imaging observations of the symbotic long-period Mira variable 
R Aquarii (R Aqr) at near-infrared and radio wavelengths. The near-infrared 
observations were made with the IOTA imaging interferometer in three narrow- 
band filters centered at 1.51, 1.64, and 1.78 /an, which sample mainly water, 
continuum, and water features, respectively. Our near-infrared fringe visibility 
and closure phase data are analyzed using three models, (a) A uniform disk model 
with wavelength-dependent sizes fails to fit the visibility data, and is inconsistent 
with the closure phase data, (b) A three-component model, comprising a Mira 
star, water shell, and an off-axis point source, provide a good fit to all data, (c) 
A model generated by a constrained image reconstruction analysis provides more 
insight, suggesting that the water shell is highly non-uniform, i.e., clumpy. The 
VLBA observations of SiO masers in the outer molecular envelope show evidence 
of turbulence, with jet-like features containing velocity gradients. 

Subject headings: stars: AGB, stars: asymmetric stars, stars: circumstellar 
shell, stars: surface features, stars: non-radial pulsation, stars: mass loss, stars: 
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individual (R Aquarii, R Aqr), binary: symbiotic, technique: long baseline inter- 
ferometry, technique: closure phase 



Introduction 



R Aqr is the brightest known symbiotic system (IWhitelock et al.lll983l ) consisting of a 
Mira primary and a white dwarf companion. The companion has not been directly detected 
yet because of the presence of nebular emissions and circumstellar material around R Aqr 
itself. The pulsation peri od of R Aqr is 387 day s. The light curves show reduced amplitude 
lasting for several cycles (IMattei fc James! Il979l ). 



Miras with close companion white dwarfs usually are classifi ed as symbiotic systems 
(lAllenlll984l : lwhitelockl[l987l : lLuthardtlfl992l : iBelczvhski et al.lboooh . A few Miras are known 
to have companions but are not (or are only very mildly) symbiotic s ystems; this includes o 
Cet = Mira, with a probable W D companion in a multi-century orbit (IReimers Sz Cassatella 



19851 : IWood fc Karovskall2004h . 



The occasional appearance of bright emissio n lines in the visib le spectrum signifies 
the presence of a hot companion around R Aqr (iMerrilll Il92ll . Il950l ). The observed wide 
and flat minima of the light cur ve of R Aqr was interpreted as the signature of a companion 
around this target (lMerrilllll956l ). Ult raviolet observations from IUE have bee n interpreted as 
evidence of a white dwarf comp anion (IM ichali tsianos. Hobbs fc Kafatoslll980l ) . The observed 



broad infrared excess at 12/im (IKenyon et al 



19881 ) could also be interpreted as evidence for 



the presence of dust heated by a hot companion; alternatively, dust shells with different 
particle sizes or composition could also explain the observed broad infrared excess at 12/im. 

An early suggestion for an orbital period of 27 yrs from radial velocity measurements 
(1 JMerrilll 1 1 9 5 Of ) was not supported by later observations. The observed decrease in the bright- 
ness of the visible and near-infrared light curves of R Aqr during 1975-78 and on two previous 
occ asions was interpreted in terms of an obscuring circumstellar dust clump eclipsing the tar- 
get (IWillson. Garnavich fc Matteill 19811 ) . These authors sugge sted that R Aqr is an ec lipsing 
binary system with an orbital period of 44 years. Subsequently, IWhitelock et al.l (119831 ) found 
a similar eclipsing signature on 6 symbiotic stars including R Aqr, however, the frequency of 
occurrence of these decreases in symbiotic stars (3 or 4 out of 6 well studied cases) with such 
long periods suggests that eclipsing by a secondary body probably is not the li kely mecha 



nism causing the observed long term amplitude changes of a few tens of years (IFeast et al. 



19831 ). However, if the Mira wind geometry is altered by the presence of the white dwarf, 
or by two winds colliding, then the eclipse could correspond to observing through the un- 
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modified Mira wind on the side away from the white dwarf. Alternat i vely, episodic dust 
emission observed around Mira and AGB stars (IDanchi fc Besterl Il995l ; iBester et al.l Il996l ) 
could explain these observed brightness modulations. A nova explosion on the surface of the 
white dwarf companion followed by dust condensation could also explain the observed long 
term amplitude changes. R Aqr could very well be a potential candidate of becoming a slow 
nova such as RR Tel which is a symbiotic system with a Mira primary and infrared excess. 

In this paper we use the word "asymmetry" to mean that part of the 2-dimensional 
brightness distribution which cannot be made symmetric with respect to a reflection through 
a point. Thus, for example, an elliptical uniform disc or an equal-brightness binary system 
are both symmetric, but a bina ry with unequal bri ghtness or a star with an off-centered 
bright /dark spot is asymmetric (IRagland et al.ll2006l ). 



Departures from circular symmetry are known in AGB stars from various hig 



h an- 



gular resolution observations 


(Karovska et al. 199 


1; Wi 


lson et al. 1992: 


Haniff et al. 


1992; 


Richichi et al. 1995: 


Raeland 


1996: Weisrelt et al. 


1996: 


Karovs 


ca et al. 


1997: Tuthill et al. 


1997; 


Lattanzi et al 


Il997: Wittkowski et al 


1998: 


Tuthill et al. 


1999, 


2000: Hofmann et al. 


2000; 


Thompson et al. 2002; ] 


Vlonnier et al. 


2004a; 


Weiner et al. 


2006; 


Ragland et al. 


2006). 



However, no consensus exists as to the mechanis m that would cause such departures from 
apparent circular symmetry (IRagland et al.ll2006l ). 



Observations of dust shells around AGB stars in the infrared 


(Danchi et al. 1994: Lopez et al. 


1997: Tuthill et al. \ 


2000: 1 


vlonnier. Tuthill & Danchi 200ol: W 


'einer et al. 20061) an 


i SiO 


maser shells at radio /millimeter wavelengths ( 


Diamond et al. 


1994: 


Greenhill et al. 


1995; 


Diamond & Kemball 


2003: 


Cotton et al. 2004: 


Soria-Ruiz et a. 


. 2004 


; Cotton et al. 


2006) 



also show asymmetries. Again, the connection between apparent surface features (or com- 
panions) and the morphology of the dust or SiO shells has not been established. 



Recent studies of astrophysical iets around a few AGB stars ( 


Kelloe;e. Pedeltv & Lvon 


2001; 


Imai et al. 


2002: 


Sahai et al. 


2003: 


Sokoloski & Kenvon 


20031: 


Brocksopp et al. 


2004) 



from radio, X-ray or Hubble Space Telescope (HST) observations suggest that those stars 
showing substantial asymmetry may all have a low mass stellar companion accreting mass 
from the AGB primary. 

The R Aqr system has complex nebular structures at arcsecond and ar cminute scales. In 



addition, X-ray, UV, vis i ble and radio observations show astrop hysical jets (IKellogg. Pedeltv fc Lyon 



200ll : iHollis et al.lll99ll : ISolf & Ulrichlll985l : iHollis et al.lll985[ ). These large scale structures 
have no influence on our high angular resolution multi-wavelength observations of the central 
region (size ~ 0.4 arcsec) of this complex system reported in this paper. 



Recently, we reported the results from the first phase of our program in which we 
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surveyed 56 evolved giants including 35 Mira stars , 18 SR variables an d 3 Irr variables 
looking for asymmetry in their brightness profiles ( fjRagland et al.ll2006l ). hereafter Paper 
I). We found that about 29% of the AGB stars show non-zero (or non-7r) closure phase 
signifying the asymmetric nature of the brightness profiles. R Aqr is the lone symbiotic star 
with positive asymmetry detected by our AGB survey. In this paper we report followup 
observations on R Aqr to characterize the observed asymmetry in this symbiotic system. A 
second objective is to detect a possible water shell around R Aqr and its possible role in the 
asymmetry detection. The reasons for combining our near-infrared observations with SiO 
maser observations are the following: (1) probably the water shell emissions detectable in the 
near-in frared and the SiO emissio ns detectable in the radio arise from the same circumstellar 
region (jCotton et al. 2004 . 2006 ). (2) the SiO maser observations give us some ideas about 
the distortions from spherical symmetry we might expect in the near-infrared, and also some 
velocity indications. 



Observations 



2.1. Infrared Optical Telescope Array (IOTA) Observations 



The observations reported here were carried out at the IOTA 3-telescope array equipped 
with an integrat ed-optics beam-combiner, IONIC , operating; in the H-band (1 -6 5 ittm) atmo- 
spheri c window (ITraub et al.l (120041 ) , iBerger et al.l (120041 ) , ISchloerb et al.l (120061 ) , iRagland et al 
(|2006h ) from 30 October to 16 November 2005. 



The near-infrared spectrum of R Aqr is dominated by water-vapor absorption bands 
(IWhitelock et al.l 119831 ). We used three sub-H-band filters to isolate the continuum and 
absorption-band part of the spectrum in order to probe the possible water shell around R 
Aqr. Our sub-H-band filters are H# = 1.51 /xm (water molecular band on blue side of center), 
He = 1.64 /im (near-continuum), and H# = 1.78 /zm (water molecular band on red side of 
center); the FWHM values are 0.09, 0.10 & 0.10 /im respectively. The bandwidths are a 
good compromise between sensitivity and spectral purity. The transmission curves for these 
narr owband filt e rs are shown in Figure U along with a n H band spectrum of R Aqr taken 
from lwhitelockl (b()07t ) (see also Iwhitelock et al.l (|l983h ). 
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Fig. 1. — The transmission curves of the narrowband filters used in this study are presented 
here. The three curves from left to right corresponds to the three wavelengths: = 1.51 
/jiii, H<7 = 1.64 jim &: H B = 1.78 /im. Also shown by square symbols is the H band spectrum 
of R Aqr from Iwhitelockl J2007h . 
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Typically, 4 minutes of program star observations were followed by nearby calibrator 
observations under identical instrumental configurations. On each star, we record 4 sets of 
data files each containing 200 scans. A scan consists of changing the optical path difference 
between two beams by roughly 75 /zm in saw-tooth form. We then take 400 scans of shutter 
data for calibration. The shutter data sequence consists of allowing only one beam at a 
time (telescope A, B and then C) and at the end blocking all three beams. Each scan takes 
about 100 ms. The field of view of the instrument in the sense of star light coupled in to the 
beam-combiner could be estimated as ~ A/D = 0.8 arcseconds for D=0.45m, the aperture 
size of the IOTA telescopes, and A = 1.65 xl0 _6 m, the wavelength of observations. The 
instrument field-of-view is further restricted by the fringe scan length to ~ tan -1 (A ^/ B max ) 
= 0.4 arcsec for B max = 26.2995m, the maximum projected baseline of the IOTA observations 
and A op d = 50 xl0~ 6 m, the smallest fringe scan size of the data block used for bispectral 
analysis. 

Interferometric fringes were recorded simultaneously on all three baselines of the IOTA 
enabling three visibilities and one closure phase measurements for each set of fringe packets. 
Observations were taken at three IOTA configurations namely, A35B15C10, A25B15C05 
and A25B05C05, providing better (u,v) coverage for the science target; where letters A, B 
and C identify the three telescope apertures and the integers denote the distance of the 
corresponding aperture in meters from the intersection of South-East and North-East arms. 
More details of the observations and instrumental configurations are given in Paper I. The 
pulsation period and phases of R Aqr during our observations are derived from AAVSO and 
AFOEV databases. A light curve of R Aqr around our interferometric observations is shown 
in Figure [21 
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Fig. 2. — The dates of observations from Table 1 (near JD 2443700) and Figures 141 [TBI and [141 
are shown on a visual light curve for R Aqr. The visual data were provided by the AAVSO, 
www.aavso.org. 
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2.2. 



VLBA Observations 



The observations of R Aqr were made as part of a more extended VLBA monitoring 
program of targets being imaged with IOTA and will be reported on more fully in a separate 
publication. The SiO maser observations of R Aqr were made using the NRAO Very Long 
Baseline Array (VLBA) on 16 September 2006. Data were recorded in 2 x 4 MHz wide 
channels in each right- and left-hand circular polarization centered on the J=l-0,i/=2 and 
J=l-0,^=1 SiO maser transitions at 42.820587 and 43.122027 GHz. Correlation of the data 
was performed at the NRAO VLBA correlator in Socorro, NM, USA and produced 128 
channels in all polarization combinations in each of the two transitions. Strong continuum 
sources near each star were observed to serve as delay, bandpass and polarization calibrators. 
R Aqr was observed in 5 observations of 16 minutes duration each. The quasar J2334+0736 
was used clS db delay calibrator for R Aqr. 



The recorded interferograms were reduced with an IDL software package. Details are 
given in Paper I. In brief, squared visibilities and closure phases are estimated from the 
three simultaneously recorded fringes through Fourier analysis. The measured fringe power 
is proportional to the squared visibility and the proportionality constant is calibrated using 
a nearby calibrator of known visibility. The closure phase is estimated as the phase of the 
bispectrum of the three fringes. 

The calibrator used for our observations of R Aqr is HD 787. The adopted angular size 
for this calibrator is 2.43 mas. Two more calibrators, namely, HD 4585 and HD 221745, 
were observed in order to monitor the stability of the instrument for squared visibility and 
closure phase measurements. The adopted angular sizes for these calibrators are 1.36 mas 
and 1.38 mas respectively. The reduced data are presented in Table [T] along with baseline 
information. 

Typical one-sigma formal errors in the mean values of our uncalibrated closure phase and 
V 2 measurements are ~ 0.2° and ~ 1% respectively. The formal errors are estimated from the 
scatters of the 200 fringe scans. We estimate calibration errors by observing calibrators under 
same observing condition and calibrating one calibrator with the other after accounting for 
the finite sizes of both calibrators. The estimated one-sigma calibration error is ~ 4% for 
the V 2 measurements and ~ 0.8° for closure phase measurements (Figure [3]). In addition, 



3. 



Data Reduction 



3.1. 



IOTA Data Reduction 
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there is another source of error in the V 2 measurements while subtracting background power 
from the total power in order to estimate the fringe power. Typical value for this additive 
error is ~ 2%. Further discussi o n on our squared visibility and closure phase measurements 
can be found in iRagland et al.l (120041 ) & Paper I. 
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Fig. 3. — Shows the calibrated squared visibility and closure phase measurements of several 
calibrator-pairs observed during 30 Oct. - 18 Nov., 2003, under identical observing configu- 
ration in order to estimate calibration errors. The estimated calibration error in the squared 
visibility measurements is ~ 4% and that of the closure phase measurements is ~ 0.8°. 



Table 1. Log of R Aqr observations 



Baseline BC Baseline AC Baseline AB 



Date 


4> 


A (fim) 


Bi (cm) 


PAi (°) 


v? 

1 


B 2 (cm) 


PA 2 (°) 


v? 

2 


B 3 (cm) 


PA Z (°) 


V? 

3 




30 Oct. 


0.06 


1.51 


1479.9 


-32.1 


0.248 ± 0.05 


1627.9 


-134.9 


0.173 ± 0.04 


2431.0 


-171.3 


0.040 ± 0.04 


-4.8±0.9 


31 Oct. 


0.07 


1.51 


1160.6 


-3.4 


0.365 ± 0.05 


2427.4 


-127.9 


0.038 ± 0.04 


3229.3 


-145.1 


0.003 ± 0.04 


-17.9±0.9 


03 Nov. 


0.07 


1.51 


1317.3 


-24.5 


0.310 ± 0.05 


2014.5 


-128.8 


0.080 ± 0.04 


2665.5 


-157.4 


0.019 ± 0.04 


-10.9±0.9 


04 Nov. 


0.08 


1.51 


1313.4 


-24.2 


0.330 ± 0.04 


2023.5 


-128.7 


0.078 ± 0.04 


2673.5 


-157.1 


0.017 ± 0.04 


-11.4±0.8 


10 Nov. 


0.09 


1.51 


1473.7 


-44.8 


0.228 ± 0.04 


1216.8 


-137.5 


0.333 ± 0.04 


1953.8 


173.6 


0.089 ± 0.04 


-0.9±1.1 


12 Nov. 


0.10 


1.51 


1006.3 


-30.2 


0.428 ± 0.05 


1904.2 


-127.6 


0.077 ± 0.04 


2264.8 


-153.7 


0.035 ± 0.04 


-7.8±0.9 


13 Nov. 


0.10 


1.51 


1440.2 


-44.7 


0.247 ± 0.04 


1299.9 


-134.7 


0.226 ± 0.04 


1940.7 


177.3 


0.083 ± 0.04 


-2.1±0.9 


16 Nov. 


0.11 


1.51 


482.1 


2.4 


0.858 ± 0.07 


1809.1 


-127.3 


0.111 ± 0.04 


2149.2 


-137.3 


0.057 ± 0.04 


2.3±0.9 


30 Oct. 


0.06 


1.64 


1164.0 


-5.1 


0.597 ± 0.07 


2409.3 


-127.7 


0.090 ± 0.04 


3190.8 


-145.6 


0.014 ± 0.04 


-10.0±0.9 


31 Oct. 


0.07 


1.64 


1479.9 


-32.1 


0.357 ± 0.05 


1627.9 


-134.9 


0.345 ± 0.04 


2431.1 


-171.3 


0.063 ± 0.04 


-3.8±0.9 


04 Nov. 


0.08 


1.64 


1523.0 


-33.4 


0.348 ± 0.05 


1519.9 


-137.8 


0.347 ± 0.04 


2404.4 


-175.7 


0.097 ± 0.04 


0.6±0.9 


06 Nov. 


0.08 


1.64 


1315.6 


-24.4 


0.369 ± 0.06 


2018.5 


-128.7 


0.163 ± 0.04 


2669.0 


-157.3 


0.038 ± 0.04 


-9.2±0.8 


09 Nov. 


0.09 


1.64 


1170.7 


-38.8 


0.470 ± 0.05 


1744.2 


-127.5 


0.296 ± 0.04 


2078.3 


-161.8 


0.138 ± 0.04 


-2.9±0.9 


10 Nov. 


0.09 


1.64 


1316.6 


-43.0 


0.423 ± 0.05 


1541.9 


-129.6 


0.319 ± 0.04 


1966.5 


-171.5 


0.171 ± 0.04 


-0.4±0.9 


12 Nov. 


0.10 


1.64 


1307.2 


-42.8 


0.423 ± 0.04 


1557.0 


-129.3 


0.321 ± 0.04 


1971.9 


-170.8 


0.151 ± 0.04 


-3.5±0.9 


15 Nov. 


0.11 


1.64 


499.5 


-13.2 


0.834 ± 0.05 


1507.3 


-130.1 


0.336 ± 0.04 


1789.4 


-144.5 


0.219 ± 0.04 


0.5±0.8 


30 Oct. 


0.06 


1.78 


1289.1 


-22.5 


0.453 ± 0.06 


2079.6 


-128.2 


0.139 ± 0.04 


2727.4 


-155.3 


0.037 ± 0.04 


-11.7±0.8 


31 Oct. 


0.07 


1.78 


1270.7 


-21.0 


0.408 ± 0.05 


2122.2 


-127.9 


0.110 ± 0.04 


2771.8 


-153.9 


0.050 ± 0.04 


-12.2±0.9 


03 Nov. 


0.07 


1.78 


1484.8 


-32.3 


0.344 ± 0.05 


1615.9 


-135.2 


0.309 ± 0.04 


2427.4 


-171.8 


0.086 ± 0.04 


-1.0±0.8 


04 Nov. 


0.08 


1.78 


1167.6 


-6.5 


0.483 ± 0.07 


2393.6 


-127.6 


0.082 ± 0.04 


3159.0 


-146.0 


0.015 ± 0.04 


-7.1±1.1 


04 Nov. 


0.08 


1.78 


1161.7 


-4.1 


0.437 ± 0.05 


2420.9 


-127.8 


0.073 ± 0.04 


3215.0 


-145.3 


0.012 ± 0.04 


-7.1±1.0 


06 Nov. 


0.08 


1.78 


1169.9 


-7.2 


0.448 ± 0.06 


2385.1 


-127.5 


0.081 ± 0.04 


3143.0 


-146.3 


0.014 ± 0.04 


-9.3±0.9 


10 Nov. 


0.09 


1.78 


1069.8 


-34.2 


0.775 ± 0.09 


1849.6 


-127.3 


0.199 ± 0.04 


2185.9 


-156.6 


0.117 ± 0.04 


-3.6±0.8 


15 Nov. 


0.11 


1.78 


483.9 


4.9 


0.821 ± 0.07 


1844.9 


-127.3 


0.196 ± 0.04 


2199.4 


-136.7 


0.091 ± 0.04 


-3.5±0.8 
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3.2. VLBA Data Reduction 



The VLBA SiO maser data were reduced as described in ICotton et al.l (12004 ) and 



Cotton et al.l (120061 ) using a combination of the AIPS and Obit software packages. Ini- 
tial amplitude calibration used measured system temperatures, weather data and standard 
antenna gain curves. Group delay measurements from the continuum calibrator were ap- 
plied to the stellar data and a spectroscopic data channel in each transition was used to self 
calibrate the phases. Instrumental polarization was determined from the weakly polarized 
continuum calibrators and the linear polarization angle from comparison with calibration 
measurements made with the VLA. Using the phase calibration derived from a single chan- 
nel, it is then possible to image each spectral channel and the channel images will be aligned 
with each other. Thus, internally consistent images can be derived but the location relative 
to the star (not visible in the images) is not measured. The restricted UV-coverage of these 
snapshot images coupled with the large angular extent of the masers required the constrained 
CLEAN available using the Obit package task Imager. The resultant images are shown with 
velocity coded as color in Figure HI The relative location of the star and the masers was 
estimated assuming that the innermost masers lie on a circle or ellipse centered on the star. 
The resolution of the derived images is represented by an elliptical Gaussian of 550 x 160 
fisec with the major axis at position angle -12°. 
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Fig. 4. — Left: SiO J=l-0,v =2 maser in R Aqr on 16 September 2006. Color indicates 
velocity, red=-16.0 km/s, blue=-26.7 km/s, and the displayed intensity is proportional to 
the square root of measured intensity. The yellow circle has the diameter of the H 2 shell 
from Table 2 (26.7 mas) and is concentric with the maser ring which is assumed to be 
approximately centered on the star. 
Right: As Left but for J=l-0,u =1 transition. 
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4. Discussion 

The reduced infrared data were fitted with various models and the derived parameters 
are discussed in this section. The parameters of these models are presented in Table [2j We 
start with a very simple model and progress towards more complex ones. Our new non-zero 
closure phase measurements presented in this paper are consistent with our earlier finding 
(Paper I). We will use this information in the more complex models. 



4.1. Wavelength Dependent Uniform Disk Sizes 

The squared visibility data taken at three wavelengths are fitted separately with a uni- 
form disk (UD) model and the wavelength dependent sizes are derived. We knew from our 
earlier work (Paper I) that R Aqr has an asymmetric brightness distribution and hence the 
UD model is not the best one for this object. Nevertheless we started this way in order to (1) 
compare our results with the UD values available in the literature, (2) investigate the role 
of possible water shell around R Aqr through wavelength dependent sizes. Moreover, the 
strength of detected asymmetric features are typically a few % of the strength of the sym- 
metric component (Paper I), implying that the asymmetric features would only mariginally 
bias the uniform disk size estimation. 

The derived wavelength dependent angular sizes are shown in Figure [5] (Top Row) . The 
derived angular diameters are 15.4±0.2 mas, 13.8±0.2 mas and 14.8±0.2 mas in the 1.51, 
1.64 and 1.78 fim bands respectively. Thus, the size of R Aqr is 12% and 7% larger at 1.51 
fj,m and 1.78/im repectively than at 1.64/iin. The reduced chi-square value for this model fit 
to the squared visibility data is 2.6. 

The near-infrared spectrum of R Aqr shows strong absorp tion bands due t o wat er 



molecules in the wings of the near- infrared bands (Figured] als dWhitelock et aJj (119831 )). 
Hence, we attribute the apparent wavelength dependent UD sizes to the presence of a water 
shell around the Mira star. 
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Fig. 5. — Top Row: Squared visibility data points are shown as diamonds along with 
the error bars. UD models fitted to the data points are shown in dotted lines. Middle 
Row: Squared visibility data simultaneously fitted with a wavelength-independent UD star 
surrounded by a water molecular shell. While optical depths at 1.51 and 1.78 /im are treated 
as free parameters, the optical depth at 1.64 /im is fixed at zero. Bottom Row: Same as 
middle row except that optical depths at all three wavelengths are treated as free parameters. 
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4.2. Wavelength Independent UD stellar size and Water Shell 



We use a two-component model, a wavelength-inde pendent UP star surrounded by a 
water shell whose optical depth is wavelength dependent (IScholal200ll ). The model visibility- 
square and closure phase are estimated through Fourier- Transform of the projected one- 
dimensional intensity profile. For simplicity, the water shell is assumed to have negligible 
radial thickness. The sizes of the star and the water shell, and the wavelength dependent 
optical depths of the water shell are derived through a model fit to the observed multi- 
wavelength squared visibility data. Assuming that the optical depth of water molecules at 
1.64/im, the pseudo-continuum band, is negligibly small, we fixed this parameter to zero in 
the model fit. The results are shown in Figure (Middle Row). Similar interpretation of 



a wat er shell present around t he centr a 



stars (IMennesson et al.l (12002 ): IWeinerl (120041 ): IPerrin et all (120041 ) ). 



star can be found in the literature for other Mira 



The free parameters of the model fit are wavelength-independent stellar angular diam- 
eter, angular size of the shell and optical depths at 1.51 and 1.78/mi. The optical depths 
at 1.64/im is assumed as zero. Stellar effective temperature is constrained from the appar- 
ent bolometric magnitude, mbol = 2.23 with an amplitude of variation of 0.91 magnitudes 
Jwhitelock. Marang fc Feastlhoooh . and the angular diameter using 



Log(T t 



eff) 



4.22 - 0.1m bol - 0.Blog((f)) 



(IRidgway et al.lll980l ). The shell radiative equlibrium temperature is constrained, using 



rp4 

shell 



2 eff 



■'star 



->sh 



'star 



Kpdr 



( 1 Chandr asekharl (119341 ): iBowenl (119881 )). Geometrically thin shell approximation would give, 



rp4 

shell 



rp4 




T~shell 



where r s h e u is the mean opacity for the shell across the part of the stellar spectrum containing 
most of stellar luminosity. 

The derived wavelength-independent stellar angular diameter is 13.2±0.1 mas and the 
derived shell diameter is ~ 2.430*. This yields an effective temperature of 2737K with 
an uncertainty due to the bolometric variability of about 11%. The radiative equilibrium 
temperature at 2.43 stellar radii, the derived shell position, is 1568K assuming T s h e ii=0. 5. 
This is a reasonable temperature for a high abundance of molecules of CO, CO2, and H 2 
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( ITsujil 119731 ). The derived optical depths of the water shell are 0.13 and 0.07 at 1.51/im 
and 1.78/zm respectively, and the average across the stellar flux should be somewhat higher 
since we do not observe where these molecules have their peak opacity (i.e. in between our 
atmospheric windows). The mean optical depth is very unlikely to be below 0.1 or 0.2, 
corresponding to shell temperatures of 1337 - 1406K. It is also unlikely to be 1.0 or higher, 
as then the shell temperature would be > 176 IK and the abundances of the molecules would 
then be lower. In fact there is probably a self-regulating mechanism operating to keep the 
optical depth below 1, since higher values lead to higher temperatures that decrease the 
abundances of the species providing the opacity. 



4.3. Importance of closure phase measurements 

So far, we have used only squared visibility data in our model fits. In this section, we 
discuss how closure phase measurements could constrain our model fits. 

The closure phase measurements of stars (with a fairly sharp edge such as a uniform 
disk) have two kinds of high angular resolution information: (1) A non-zero or non-7r closure 
phase signifies an asymmetric brightness profile (Paper I). Asymmetric here means structure 
which is not rotationally symmetric. (2) The rough magnitude of closure phase (whether 
close to zero for a small star or tt for a large star) constrains the size of the central star 
and thereby could require additional extended circumstellar components to fit the squared 
visibility data. 

Centrosymmetric targets should give a closure phase of eit her zero or ± tt dep ending 



on how many baselines are beyond the first, second, etc., nulls (IRagland et al.l 120061 ) . The 
measured visibilities as a function of baseline and the closure phase data suggest that all 
measurements reported in this article are taken well within the first null of R Aqr. Hence, 
the longest baseline data (see Tabled]), taken on 31 Oct 2005 at 1.51/im with B maa .=32.293 
m that gives a closure phase of 17.9° could provide an upper limit on the possible uniform 
disk size for the central star. In order for the closure phase measurement to be close to zero 
rather than tt, the diameter of the central star should be less than 1.22 A/B, i.e. the diameter 
must be < 11.8 mas. Thus, the closure phase data clearly shows that we have overestimated 
the sizes in the previous two models by at least 13% with the H 2 shell model (Section 4.2) 
and in the range 17% to 31% with UD model (Section 4.1). 

Thus, the potential of closure phase measurements is not only in detecting asymmetry 
but also in detecting extended symmetric features around the central star. This also suggests 
that the UD sizes of AGB stars derived from limited squared visibility measurements should 
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be treated with caution. 



4.4. Improved two component model: Star surrounded by a water molecular 

shell 

In Figure (bottom row), we present a two component model of a central star with 
wavelength independent size surrounded by a thin shell due to water molecules. The sizes 
of these two components and the wavelength dependent optical depth are derived through 
non-linear least-square fit to the observed multi-wavelength squared visibility data. The only 
difference from the earlier model (Section 4.2) is that the optical depth at 1.64/xm is also 
treated as a free parameter. Still, we have not used our closure phase data, but the derived 
size for the central star is now consistent with our upper limit on the stellar diameter from 
closure phase values (see Section 4.3). 

The derived wavelength independent stellar angular diameter is 11.9±0.2 mas and the 
derived shell diameter is 2.320*. As before, the stellar effective temperature is derived as 
2881 K and that of the shell is estimated as 1691 K assuming r s h e u=0.5. The derived optical 
depths of the shell are 0.21, 0.07 and 0.13 at 1.51/xm, 1.64 /zm and 1.78 /im respectively. The 
derived stellar size is consistent with the closure phase measurements. We estimate relatively 
large optical depths in all three wavelengths. The non-zero optical depth at 1.64//m could 
be attributed to a limb darkening effect, scattering of molecules and dust around the star 
and possible absorption lines in the 1.64/im filter band. The model over fits the squared 
visibility data with a reduced chi-square of 0.8. 

4.5. Three component model: Star, Water Shell and an off-axis compact 

feature 

The next step in our analysis is to account for the non-zero closure phase measurements. 
As shown in Paper I, the non-zero closure phase could only be explained with an off-center 
additional component. It could be an off-center (bright or dark) spot on the stellar surface 
or a companion or dust clump. The simplest possible model would be to add a point source 
to the existing two-component model. 

Figure [6] shows the results of this three component model. Addition of this third com- 
ponent marginally made changes to the star and shell parameters. The parameters of this 
three component model are given in Table [2] (Model 4). This three-component model si- 
multaneously fits the squared visibility and closure phase data with a reduced chi-square of 
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1.15. We also found an alternate solution through Monte Carlo simulation by systematically 
changing the position of the unresolved compact feature; this model is given in Table [2] 
(Model 5). The derived optical depths of the shell are 0.23, 0.08 and 0.13 at 1.51/im, 1.64 
jum and 1.78 /im respectively for Model 4, and 0.26, 0.09 and 0.14 at 1.51/mi, 1.64 /zm and 
1.78 /im respectively for Model 5. The primary difference between the two models (Model 
4 & Model 5) is the position of the off-axis model feature. These de rived optical depth s are 



consistent with that reported in the literature for other Mira stars (IPerrin et all 120041 ) . 



The derived wavelength independent stellar angular diameter is 11.2 mas and the derived 
shell diameter is 2.380*. The position angle of the fitted model feature (off-axis point source) 
is 131° which is roughly orthogonal to the axis of the astrophysical jets detected in the x-ray, 
UV, visible and radio wavelengths and the angular separation of the off-axis compact feature 
with respect to the central star is 11 mas. The position angle and angular separation of the 
off-axis compact feature for the alternate model (Model 5) are 266° and 7.5 mas respectively. 



Tuthill et all (120001 ) observed R Aqr at near-infrared wavelengths using aperture mask- 
ing interferometry at the Keck-I telescope with an angular resolution of ~ 15 mas. These 
authors did not detect any companion with a relative brightness of < 5 magnitudes. This is 
probably consistent with any of the models in Table 5, since a 2-5% feature at a separation 
of 0.5-0.7 times the angular resolution is not likely to have been detected. 
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Fig. 6. — Squared visibility and closure phase data fitted with a three-component model: 
a wavelength independent UD star surrounded by a molecular shell and a off-axis model 
feature is shown here. The data points along with the error bars are shown in cross symbols 
and the model values are shown as triangles. 



Table 2. Results of the four model fits. 
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1 These values are not used for chi-square optimization (see text). 
4>\ and <j>sh ar e the angular diameters of the central star and the molecular shell respectively. 50 is the angular separation between the central star and off-axis model 
feature. 
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4.6. Accretion Disk Hypothesis 

In this section we focus on whether or not the fitted model feature in the IOTA observa- 
tions can be accounte d for by the white dwarf companion that has been previously discussed 



in the literature (e.g. jHollis et al.l (119971 )) or by the effect of Roche lobe hot spot. First, we 
discuss the luminosity ratio between the white dwarf companion and the AGB star. Next, we 
calculate the effect of the accretion disk to determine its luminosity contribution (4.6.2), and 
a number of other relevant physical effects including mass loss on the accretion luminosity 
(4.6.3), heating of the AGB star from the white dwarf (4.6.5), heating of the accretion disk 
from the AGB star (4.6.6), and the luminosity of the hot spot of the accretion disk (4.6.7). 
Understanding the magnitude of the infrared emission due to these physical effects allows 
us to determine whether or not the companion feature is consistent with a white dwarf star 
with the orbital parameters in the literature, or if it could be an additional white dwarf 
companion, or if the companion feature is really due to a physical effect of the AGB star 
itself, such as a hot spot near the edge of the star. 



4-6.1. White dwarf Luminosity 



We begin this discussion with a simple calculation of the luminosity of the white dwarf 
companion compared to that of the red giant based on parameters for these stars in the 
literature. We assume the red giant star ha s a radius of 300 R , a temperature of 2800 



K, giving it a total luminosity of 5 x 10 3 L p (IKafatos et al.l Il986l ; iHollis fc Koupelis 1 12000 



"1 • 

Bureare 


la et al. 


Il9 




Makinen et al. 


2004) 


work of 


Hollis et al. 


(1997h- 



Two different sets of parameters for the companion have been used in the literature. One 
choice is that the companion has a temperature of around 40,000 K and a diameter of 0.1 
R , which would be classified as a subdwarf star, with a total luminosity of approximately 
10 L . The other choice is a white dwarf star, with a temperature of about 130,000 K, and 
a diameter of 0.014 R , also giving approximately 10 L Q for the total luminosity. 

We begin by calculating the ratio of the luminosity at 1.65 fim of the companion to 
that of the red giant star. The result is that the ratio for the subdwarf case is about 10~ 5 , 
whereas for the white dwarf the ratio is 7.1 xl0~ 7 , which is substantially smaller. Given that 
the ratio of luminosities for the observed asymmetry corresponds to about 2% of the total 
1.65 fim luminosity of the red giant star, it is clear that the asymmetry cannot be due to 
the detection of the companion itself. Hence the accretion disk around the white dwarf star 
or subdwarf companion must contribute the bulk of the luminosity at this wavelength to 
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explain the observational results. For the remainder of this section we explore a number of 
physical effects associated with the companion and its accretion disk that could explain the 
observations. 



4-6.2. Accretion Luminosity 



The ac cretion luminos ity can be readily calculated based on the results of lLynden-Bell &: Pringle 



(1974 


) and 


Prinele 


(1981) 



For an optically thick, physically thin, flat disk, the luminosity at the earth is given by; 



f*° ut 2n RnB„{T(R)]dR 



D 2 



where 



T(R) 



( 3GM*M 
I 8ttct 



1/2' 



1/4 



(2) 



is the temperature in the accretion disk (jPringle Ill98ll ) and D is the distance to the Earth. 
The quantity B u is the Planck function, M* is the mass of the white dwarf, R* is its radius, 
M is the accretion rate onto the disk, and o and G are the Stefan-Boltzmann constant and 
the gravitational constant, respectively. The maximum temperature in the accretion disk 
occurs at a radius R ma x = 1-36 R*, and is given by 



T 



0.488 



3GM*M 
%-ko-RI 



1/4 



(3) 



as discussed in ILynden-Bell &: Pringld ( 119741 ) , iPringle I (Il98ll ) , and lHartmann I (119981 ) . In 
this treatment we neglect the emission from the boundary layer, which occurs over a very 
small angular region around the surface of the white dwarf star and emits at a very high 
temperature, which would have a negligible contribution to the luminosity in the infrared 
region near 1.65 ^m. 

With the parameters for the white dwarf star and the AGB star as discussed previously 
in this section, we can calculate the ratio of the AGB luminosity to that of the accretion 
disk. The red giant luminosity is given by: 
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L AGB (v) = ir\^j B U (T*) (4) 

where R*, T* are the radius and temperature of the star, and D is the distance to the Earth. 

From these equations we calculate the ratio of the accretion luminosity at 1.65 fim to 
that of the AGB star, assuming different values for the mass loss rate from the AGB star 
and an adopted efficiency of capture of material onto the accretion disk of 10%. Thus, if the 
mass loss rate were 10~ 5 M Q yr _1 , then the accretion rate would be 10 -6 M yr _1 . Table 
[3] presents our results for varying mass loss rates. We find that the maximum contribution 
to the 1.65 /im flux from the accretion disk occurs at the highest conceivable mass loss rate 
from the AGB star, namely at 10~ 5 M Q yr -1 , and the ratio is not strongly dependent on 
the assumptions about the white dwarf star. The ratio of 1.6% (Table 3) is reasonably 
close to that of the closure phase measurements for the expected intensity ratio, so we must 
determine whether the mass loss rates are very high, as required to explain the measurements, 
or whether they are relatively low. 
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Table 3. Accretion luminosities of White Dwarf companion to AGB star at 1.65 fim. 



WD WD AGB Mass WD T a cc,max Ij AGB /^Accretion 

Radius Temperature Loss Rate Ace. rate (K) 
(Rsun) (K) (M Sun yr- 1 ) Rate 1 



0.1 


4.0x 


10 4 


10~ 5 


10~ 6 


4.8x 10 4 


1.5x 


IO" 2 








10~ 6 


io- 7 


2.3x 10 4 


3.0x 


IO" 3 








io- 7 


io- 8 


1.3x 10 4 


5.7x 


io~ 4 


0.014 


1.3 x 


10 5 


10~ 5 


IO" 6 


1.8 x 10 5 


1.6x 


io- 2 








IO" 6 


io- 7 


1.0 x 10 s 


3.4x 


10~ 3 








io- 7 


io- 8 


5.6x 10 4 


7.1x 


io~ 4 



WD accretion rate is assumed to be 10% of AGB Mass Loss Rate. 
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4-6.3. Mass loss estimates 
A variety of mass loss rates hav e been estimat e d for R Aqr. The largest rate is 8.9 x 10 



-7 



M yr 1 (Danchi et al. 1994), while ISpergel et al.l (119831 ) and Hollis et al. (1985) have much 



lower rates, 1.9xl0 -8 M yr -1 and 4.2xl0 -8 M Q yr -1 , respectively. These measurements 
were made using a variety of techniques, for example Danchi et al. (1994) used infrared 
interferometric measurements of the distribution within a few stellar radii of the star itself, 



while ISpergel et al.l (119831 ) used measurements of line profiles of CO, to infer mass loss rates. 



The largest of these gives a luminosity ratio of 0.3%, which is too small to account for 
the observed closure phases. To obtain accretion luminosities large enough to account for the 
observations, a much larger fraction of the AGB wind (than 10%) must have been captured 
by the white dwarf star than was used in the previous estimate. Alternatively, the more 
recent measurements were taken during a period of substantially larger mass loss than were 
used in our estimates in the previous section. 

Both alternatives are possible, as Danchi et al. (1994) and Bester et al. (1996) have 
shown that mass loss can be episodic and that there can be periods of increased mass loss rela- 
tive to the mean mass loss rate, that would be obtained from CO line profile measurements, 
which probe a much larger volume around a star than the interferometric measurements. 
Thus, the CO measurements probe a longer history (of the order of 1000 years) and hence 
average over periods when the mass loss could have been relatively low. 



4-6.4- Orbital parameters 

It is also possible that the fraction of the AGB wind captured could vary substantially in 
time, and this is due to the fact that the orbit is eccentric. Figure [7] displays the orbit of the 
white dwarf companion to the AGB star over the course of a complete orbital period, P=44 
years, beginning in 1974, the date of periastron, at T=2442100.0 JD. The parameters used 
in the calculation are from Hollis, Pedelty, and Lyon (1997) and they are: a, the semi-major 
axis, between 16 and 18 AU for a total system mass between 2.5 and 3 M Q ; the eccentricity, 
e=0.8; the inclination, i ~ 70 degrees (in their paper); the longitude of the ascending node, Q 
~ 90 degrees; and the longitude of periastron, u ~ ± 90 degrees. Thes e orbital par a meter s 



do not produce an orbit that matched the position angle reported in iHollis et al.l (119971 ) . 



Instead we find that we can match the reported position angle if we use u = -90 d egree s and 



= -75 degrees. Given the large uncertainty in the parameters in IHollis et al.l (119971 ). the 
adopted parameters are reasonable and consistent with all the data discussed herein (In the 
absence of radial velocity measurements, the inclination angle is subject to an ambiguity. 
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Accordingly, the inclination, i ~ 110 degrees; the longitude of the ascending node, Q ~ +105 
degrees; and the longitude of periastron, to ~ +90 degrees is an alternate orbital solution 
for the R Aqr binary system). The plot (Figure [7]) displays the orbit (solid line), location 
at four year intervals as labeled, and the AGB star (gray filled circle at the origin) and the 
white dwarf accretion disk (outline in black) foreshortened due to the high inclination of 
the orbit, assuming the accretion disk and the stars orbit in the same plane. The AGB star 
and accretion disk are also to scale, where we have adopted a distance of 200 pc. It's worth 
mentioning that the orbital parameters derived by Hollis et al. 1997 employed a number 
of assumptions based on largely circumstantial evidence and a single measurement of the 
offsets between the centroids of the SiO masers surrounding the Mira and the continuum HII 
region thought to be surrounding the white dwarf companion. While the centroids of these 
features may not be coincident with the star responsible for them, there is a good case made 
that in 1996, the companion white dwarf appeared several 10 's of mas north of the Mira. 
Unless the orbital period is much shorter than the assumed 44 years, the object detected by 
Hollis et al. 1997 is unlikely to be the fitted model feature (Model 4 & Model 5). 



Mcintosh &: Rustanl (120071 ) analyzed radial velocity data obtained from 1946 until 2007 



from a SiO maser emission, near-infrared and visual spectral lines. They determined an 
orbital period of 34.6 ± 1.2 yr, an eccentricity of 0.52 ± 0.08, and a projected semi-major 
axis of 3.5 ± 0.4 AU. Their results indicated periastron passage occurred in 1980.3 ± 0.8 
and the angle between the line of nodes and periastron is 110.7 ± 18.4 degrees. These 
parameter values would indicate Roche lobe overflow at periastron, however, these results 
are suspect as they produce an unphysical system mass of 0.043 M , which disagrees with 
the ex pected system mass of 2.5 to 3.0 M as discussed previously. iMcIntosh fc Rustan 
( 120071 ) can obtain a somewhat higher system mass if they use a much smaller inclination 
angle of 20 degrees, rather than ~ 70 degrees, as expected from other analyses. Continued 
monitoring and analyses of imaging and spectral data is needed to produce more precise 
estimates of the orbital elements for this system. 



4-6.5. Heating of primary from secondary 

There are a number of other effects that are worth exploring, given knowledge of the 
orbit of the two stars. One type of effect that has been noted in the literature on cataclysmic 
variable stars is a warming effect on the primary star due to radiation from the secondary 
star. For example, it is possible that the white dwarf star, which is illuminating the AGB 
star, could potentially heat up one side of the star, and create an asymmetry, which could 
be observable in the closure phases at 1.65 /im. Based on the work of Warner (1995), it is 
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easy to show that the fractional temperature change of the AGB star, AT/T, is: 



AT 7T (R 2 




2 



T 8 V d 



(5) 



where R 2 and L 2 are the radius and luminosity of the AGB star, respectively, and L x is 
the luminosity of the white dwarf star, and d is the distance between them. At periastron, 
Eqn. (5) gives a fractional temperature change of 2xl0~ 5 , which is a temperature change of 
0.056 K, which is a very small effect indeed. At the distance in 2004, d=26 AU, this effect 
is completely negligible. 



However, as the white dwarf star heats the AGB star, the AGB star likewise heats the 
white dwarf star itself and its accretion disk. Given the large luminosity of the AGB star, ~ 
5000 L Q , the effect of the AGB star on the accretion disk, particularly its cooler, outer parts, 
could be substantial. To understand the size of this effect, we make a simple calculation of 
a black dust grain in the presence of this luminosity: 



where R s tar and T star are the radius and temperature of the star and d is the distance from 
the star to the dust grain. We assume the AGB star has a temperature of 2800 K and a 
radius of 300 R . For a grain at 5.3 AU, we obtain a temperature of 1016 K, while at larger 
distances, such as 17 AU, the temperature is substantially lower, about 567 K, and at the 
location of the white dwarf in 2004, 26 AU, the temperature is 458 K. 

We can estimate the increase in luminosity of the accretion disk due to the warming 
effect from the AGB star. Using eqn. (1) and assuming a constant temperature from about 
0.5 of the outer radius of the accretion disk (about 1 AU), to the outer radius, we find the 
luminosity at 1.65 /im is 20 Jy or 1.3% at 5.3 AU, and 0.6 mJy or 4xl0~ 5 % at a distance 
of 26 AU. For the current epoch, the contribution to the 1.65 /xm luminosity is negligible, 
but there will clearly be a strong effect as the stars approach periastron, where the distance 
is a(l-e)=3.4 AU, where we have used the orbital parameters for the lower total mass. In 
this case, the effect is of the order of 113 Jy or 7% of the 1.65 jum flux density of the AGB 
star, so this should be observable with current and future long baseline interferometers in 
the next few years. 



4- 6. 6. Heating of secondary from primary 




(6) 
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4-6.7. Effect of Roche lobe hot spot 

Finally, another effect that could potentially be of importance observationally as the two 
stars approach periastron is the emission from a bright spot caused by Roche lobe overflow. 
A simple estimate of the flux produced by the bright spot is given in Warner (1995), pp. 
81-82. The flux, which we write as, F bs (i/), is 

where AR S the radial size of the spot, A8 is the radial angle along the accretion disk taken 
up by the spot, is the outer radius of the accretion disk. The symbols d and D are the 
distance between the two stars and the distance to the observer, respectively. At periastron, 
we find the hot spot is expected to contribute as much as 110 Jy (or 7%) at 1.65 //m, assuming 
a typical hot spot temperature of ~ 10000 K. 

We conclude that a hotspot from Roche lobe overflow is energetically capable of pro- 
ducing a feature that could, then, explain the observed closure phases. However, Roche lobe 
overflow is not expected for the orbit shown in Figure [TJ as the separation is much too large. 
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Fig. 7.— Orbit of the white dwarf companion to R Aqr based on the orbital elements 
published by Hollis, Pedelty, and Lyon (1997) with the longitude of the ascending node 
and periastron modified as discussed in Section 4.6.4, computed for the epoch of the IOTA 
observations. The Mira star, H 2 shell and the fitted model feature (at the two possible 
locations from Table 2) are also shown here. The SiO maser emission in R Aqr on 16 
September 2006 is superimposed on to this figure where color indicates velocity, red=-16.0 
km/s, blue=-26.7 km/s. 
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4.7. Clumpy water shell 

So far we have assumed that the central star and the water shell are symmetric in 
nature and attributed the observed asymmetric feature to an off-axis point source (accretion 
disk around white-dwarf companion or a bright hot spot). In this section, we explore the 
possibility of an asymmetric water shell (without a need for a close companion or hot spot) 
in explaining the observed asymmetry. We replaced the model fitted off-axis point sources 
(section 4.5) with Gaussian-shaped marginally-resolved features and found no significant 
degradation in the reduced chi-square. This signifies that a clumpy molecular shell would 
probably fit the measured data. In the absence of any a-priori knowledge on the shape of the 
asymmetric shell, our obvious choice was to go for a model independent image reconstruction 
algorithm to explore this possibility further. The limited uv coverage of our observations 
(Figure [TO]) would prevent reconstruction of the entire instrument field of view of ~ 0.4 
arcsecond. However, our Monte Carlo simulations looking for an unresolved feature around 
R Aqr (Section 4.5) did not detect any feature in the annular region with an inner radius 
of ~ 20 mas and the outer radius of ~ 0.4 arcsecond. Hence, we performed the image 
reconstruction of the central ~ 20 mas area around the central star. 



4-7.1. Direct minimization algorithm for image reconstruction 

The van Cittert Zernike Theorem states that the spatial coherence function is the Fourier 
transform of the brightness distribution of the astronomical source under certain conditions. 
Thus, in principle, an inverse Fourier transform of spatial coherence function should provide 
the brightness distribution of the target. Long baseline interferometeters measure only fringe 
visibility amplitudes (but not fringe phases) as a function of baselines and closure phases for 
all closed triangles. Hence, lack of fringe phase information challenges this simple Fourier 
inversion technique to reconstruct the images. 

Image reconstruction from measured visibility and closure phase data is a well known 



procedure in the context of speckle interferometry (I Weigeltlll99ll ). The challenges have been 
in adapting this method to optical interferometric arrays with sparse apertures. There is 
no unique way of expanding the limited baseline data to all of the bispectral domain or 
its subdomain where the unique image reconstruction of an arbitrary object is possible. 
However, the positivity constraint could be used to recover the image of relatively compact 
objects. 



The standard radio interferometry self calibration techniques assume direct measure- 
ment of the complex visibilities whereas the much shorter atmospheric coherence times of 



-32 - 



Optical/IR interferometers requires different techniques for estimating visibility squared from 
a power spectrum analysis and closure phase from a bispectrum analysis. This means that 
the usage of radio techniques are difficult to apply to our data. Therefore we developed an 
algorithm based on direct minimization procedure where the model image plane intensity is 
varied to reduce the function 



a 



- \ 2 2 



under the following constraints 



l \ i / i,j,k x 



/(r) > 0, if r < r 

/(r) = 0, otherwise, (9) 

where /(r) is the object brightness distribution and = ViVjVkexp[i(pi + <pj + <pk)] 
represents bispectral elements of the reconstructed image. The Fourier component of the 
reconstructed image for 1-th baseline is written as Viexp(i<pi). The summation is performed 
for all measured visibilities Vi and bispectral elements ViVjVkexp{i(pijk) where are mea- 
sured closure phases {(pi + (pj + pk). The measured visibilities and bispectral elements are 
weighted according to their errors a V 2 and (TB ijk - The relative contribution of pure visibilities 
and phase related values is balanced by using the constant a. The parameter r defines the 
field of view of the interferometer. 

Th e minimization p rocedure enforces positivity constraints by using the multipliers 



i_p] 

method iBertsekad (jl982l ). In this minimization procedure, the constraints (Equation [9]) are 



combined with the objective function (Equation [8]) to form a modified Lagrangian function 
of the form 



Lc = F + - ^([max(0, 7i - c/(r,))] 2 - 7 J) (10) 

3 

and the minimization problem (Equations E]& W) is replaced with a sequence of unconstrained 
minimization procedures of the form, 

minimize L Ck (f k ,^) (11) 

where {ck} is a growing positive sequence of penalty parameter c&, and minimization is 
performed for all n restored points of the image inside the FOV. And, f k (r) is the real- 
valued solution that provides the minimum of the modified Lagrange function L Ck which is 
used as the first approximation for the minimization of L Ck +1 (f, 7). 
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The new penalty parameter for the next minimization step is computed as 



c k+1 = (1 + f3)c k , 



(12) 



where (3 is a constant > 0. 

For each new penalty parameter Ck+i the Lagrange vector 7*. is recalculated as 



The vector 7 fc is set at the beginning of the iterative procedure to zeros. The uncondi- 
tional minimization of L Ck can be performed by using any known minimization method. We 
adapted gradient search method because of its simplicity. 

The conditional minimization problem (Equations [S] & E]) is incorporated into a self 
calibrating iterative loop where the output of the previous minimization is used as the 
initial approximation for the next minimization step. A Gaussian brightness distribution 
minimizing Xy2 is assumed as the initial approximation for the first minimization loop. The 
resolution of the final images is constrained by the apodization of the Fourier spectrum with 
a Gaussian function of a ~ 0.7 B max . A flow diagram of this minimization algorithm is 
shown in Figure 
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Fig. 8. — A flow chart of the image reconstruction algorithm is shown in this figure. The 
dashed-line rectangle encloses the constrained-minimization procedure. The outer loop (out- 
side this dashed-line rectangle) represents the self-calibration loop. As shown in this flow 
chart, the self-calibration loop uses the most recent values of /(r,-) as the initial values for 
f k (rj), while setting {^ k } and {c k } to initial values, namely, {7 } and {c } respectively. 



4-7.2. Reconstructed images 



The reconstructed images are shown in Figure [91 The reduced chi-squares for these im- 
ages is estimated by summing the reduced chi-squares for all visibility and the closure phase 
data points with values of Xr — 1-96, 0.85 and 1.02 respectively at 1.51/im, 1.64/xm and 
1.78/mi. For clarity, we subtracted the central star contribution from these reconstructed 
images and these are also shown in this figure. These latter images clearly show an asym- 
metric shell with at least three clumpy features. The uv coverage of our observations in 
shown in Figure [TUl 
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Fig. 9. — Reconstructed near-infrared images of R Aqr are shown here. The blue, green and 
red color images (from left to right) represent 1.51, 1.64 and 1.78 /zm respectively. The top 
row shows reconstructed images and the bottom row shows primarily contributions from the 
shell since we subtracted a Gaussian function from the reconstructed images to remove the 
stellar component. North is up, east to the left. 
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Fig. 10. — The uv coverage of our observations are shown here. The symbols "plus", "star" 
and "diamonds" correspond to the observational wavelengths of 1.51/zm, 1.64/zm and 1.78;um 
respectively. 
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Two of those clumpy features are roughly seen at positions where we obtained an off- 
axis model feature for our three component models (section 4.5). Figure [TT1 shows the radial 
and azimuthal profiles of these reconstructed images. The relative surface brightness of the 
shell with respect to the peak brightness of the image at 1.51, 1.64 and 1.78 /mi is 7.8%, 
4.4% and 6.4% respectively. The shell is brighter at 1.51 /im and 1.78 /xm than at 1.64 /mi. 
This signifies the presence of water molecules in the shell. However, non-zero shell brightness 
at 1.64 /im signifies contributions from scattering as well as dust emissions at these three 
wavelengths. 
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Radius (milliarcsec) Position angle (deg) 

Fig. 11. — Left: The radial profiles of the azimuthally averaged, reconstructed images are 
plotted here. The dot-dashed, dotted and dashed lines correspond to 1.51 pm, 1.64 fim and 
1.78 /im respectively. Right: The azimuthal profiles of the radially averaged, reconstructed 
images are plotted here. The dot-dashed, dotted and dashed lines correspond to 1.51 /im, 
1.64 /im and 1.78 /im respectively. The continuous curve corresponds to the mean values. 
The azimuthal profiles clearly show at least three peaks. Two of those features correspond 
to off-axis point source solutions of the three component model presented in section 4.5. 
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4-7.3. Simple model 

The three wavelength images are composed and shown in Figure [12j The blue, green 
and red colors represent 1.51, 1.64 and 1.78 /mi respectively. The 1.51/mi image shows a 
bipolar feature. Moreover, the features at 1.64 and 1.78/xm are somewhat consistent with 
this feature seen at 1.51/im. Interestingly the axis of symmetry of this bipolar feature is 
roughly parallel to the large scale astrophysical jets of R Aqr as well as the jets seen in 
SiO maser observations (Figure H]). We propose a simple speculative model in which the 
astrophysical jets to the north-east and south-west compress the material in these directions 
giving rise to the feature in Figure O Additional measurements are essential to validate 
this simple model. If this hypothesis is true, then these observations are key steps in the 
understanding of the formation of bipolar nebulae. 
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Fig. 12. — The three wavelength images are composed and shown in this figure. The blue, 
green and red colors represent 1.51, 1.64 and 1.78 fim respectively. North is up, east to the 
left. 
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SiO Masers 



The molecular envelopes of oxygen rich Miras frequently exhibit SiO masers (IReid fe Menton 



19981 ) in the region just interior to where the molecular gas condenses into dust. The masers 
are strongest at the locations where the line of sight is tangent to the region giving rise to 
the masers, as these are the locations where the long paths through constant radial velocity 
gas needed to develop sufficient gain in the masers occur. The masers appear as partial rings 
of bright spots centered on the star. 

The IOTA observations (table 1) are taken in the early declining phase of one pulsation 
cycle, while the VLBA data were taken in the mid-rise of the next. However, the travel 
times are such that a shock emerging at some time, t=0, through the photosphere will 
arrive at the SiO / molecular shell ab out one pulsation peri o d late r. Furthermore , moni toring; 
of several maser sh ells reported in iDiamond fe Kemballl (120031 ). ICotton et al.l (120041 ). and 
Cotton et al.l (120061 ) show a complex interaction of the maser distribution with the pulsation 
cycle with variations of the maser distribution only weakly coupled to individual pulsation 
cycles. Thus, the distribution of masers appears to be accumulated over several pulsation 
cycles; we assume that the separation of the IR and SiO observations by a cycle does not 
affect our interpretation of the results. 



Interferometric observations of R Aqr in the z/ =l and v 



tions of SiO at 43 . 1 and 42.8 GHz a re reported in Boboltz. Diamond fe Kemballl (119971) , 



=2^ 



=1-0 rotational transi 



Hollis et all (I2000h . iHollis et all (120011) . , ICotton et a 



average SiO ring diameter reported by ICotton et al. 



(120041) and ICotton et al.l ([2009 1. The 



(120061 ) for R Aqr is 31.5 mas or some- 



what larger than the H 2 shell diameter of 25.2 mas given in Table 2 or shown in Figure 



Perrin et al. 


(2004 


) and 


Cotton et al. 



(120061 ) for Mira and U Ori. 



The SiO masers shown in Figure [7] appear well outside of the H 2 shell. This suggests 
that the molecular region of the circumstellar envelope is more extended than given by the 
simple model fit or Figure [91 Figures H~3l and IT41 show velocity gradients in the outer regions 
of the molecular envelope. Since these measured velocities are primarily tangential to the 
molecular envelope, they do not necessarily represent an overall expansion of the envelope 
but may well be part of a dust expulsion event as the gas giving rise to the masers condenses 
into silicate dust. 



The observations of IHollis et al.l (120011 ) and lCotton et al.l (120041 ) show that around 2001, 
and before, there appeared to be either general differential rotation of the envelope or a 
rotating molecular ring. This feature, if interpreted as rotation, would have a period of 22 
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sin(i) years where i is the inclination of the orbit. The proj ected position angle o f the ring on 



(Cotton et al. 


2004 


)■ 


Hollis et al. 



( I2000l ) interpret this systematic motion as being the result of tidal effects on the molecular 
envelope during the periastron passage of the companion past the Mira. This systematic 
motion was repo rted in observations between 1996 and 2001 but subsequent observations in 
September 2004 (ICotton et al.l 120061 ) show no evidence for systematic rotation. 



Figure H] shows substantial differences from earlier observations given in ICotton et al. 



( 120061 ). SiO masers around AGB stars generally consist of localized features forming a partial 
ring. Turbulence will reduce lengths of the line of sight through constant radial velocity gas, 
reducing the gain and hence the brightness of the maser spots. A ring-like structure is hard 
to determine from Figure HI suggesting that the molecular envelope is more turbulent than 
usual. In Miras better studied as a function of pulsation phas e, no clear relationship has 



been found betwee n pulsation phase and SiO maser structures ( iDiamond fc Kemballl 12003 



Cotton et al.ll2006l ) and clear, if partial, rings are the norm. In particular, both transitions 
show prominent jet-like features to the north-east containing velocity gradients at position 
angles between 20 and 35 degrees. The velocity of these features becomes increasingly 
far from systemic with increasing distance from the star. Close-up views of these regions 
with polarization vectors are shown in Figures [13] and [TH The orientations of the features 
seen in these figures are quite different from the position angles (80, -100 degrees) of the 
2001 features. The polarization vectors seen in these plots are mostly perpendicular to the 
direction of the f eature indicating tha t the magnetic field is predominantly along the feature 
(see discussion in lCotton et al.l (120061 ). Thus, the range of orientations of the jet-like features 
means that it is unlikely that the magnetic fields in these features are part of a large scale, 
persistent dipolar field. 

The jet-like features seen in 2006 are approximately parallel to the large scale jet from 
the companion. In principle, the features could have been formed by the large scale jet. 
However, the 2001 jet-like features show little relationship to either the direction of the 
large scale jet or to the derived plane of the orbit of the companion shown in Figure [7J 
Furthermore, the orbit derived for the well established companion places it far from the 
AGB star just prior to 2001. Thus, any connection between the companion star or its jet 
and the SiO features is doubtful. 
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Fig. 13. — Contours of integrated Stokes' I at powers of 2 times 0.1 Jy/beam of the SiO J=l- 
0, v =2 maser in R Aqr. Superposed are vectors with length proportional to the integrated 
linear polarization and orientation of the E-vectors. The gray-scale gives the flux weighted 
average velocity whose scale is given by the wedge at the top of the plot. The resolution is 
shown in the lower left corner. 
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Fig. 14. — Like Figure [T3l except for the SiO J=l-0,z/ =1 transition. The resolution is shown 
in the lower right corner. 
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Conclusions 



In order to characterize the observed asymmetry reported in our earlier work (IRagland et al. 



20061 ). we carried out follow-up observations of selected Mira stars at near- infrared wave- 
lengths. In this article, we presented observations of the nearest known symbiotic Mira, R 
Aqr, from the IOTA imaging interferometer and from the VLBA. 

The reconstructed near-infrared images of R Aqr are reported for the first time at three 
wavelengths namely, 1.51, 1.64 and 1.78 /im. The near-infrared images suggest a clumpy 
molecular shell around the central Mira star. The observed data could also be modeled 
by a three component model consisting of a symmetric central star surrounded by a water 
shell with a radius of about 2.25R*, and an off-axis compact feature at about 2R* at a 
position angle of 131° (an alternate model exists wherein the off-axis feature is at about R* 
at a position angle of -94°). However, the comp anion feature solut ion is inconsistent with 
the orbital parameters reported in the literature (jHollis et al.lll997l ). We conclude that our 
observations are best explained with a clumpy, extended H 2 circumstellar envelope. 

The ratio of the size of water shell to the size of the central star is 2.25. The ratio of 
the size of SiO shell to that of water shell is 1.25. Thus, the masers appear at the outer edge 
of the molecular envelope as reported for other Mira stars in the literature. 

The SiO maser emission around R Aqr shows large scale turbulence in the molecular 
envelope with prominent jet-like features towards the north-east. These features are roughly 
parall el to the large scale jet from the more dista nt companion star discussed by 



( 119971 ) . However, similar features seen in 2001 (jCotton et al.ll2004l ; iHollis et al. 



lollis et al. 



200l|) have 



a very different orientation from either the jet or the inferred orbit of the companion so any 
connection between the SiO features and the companion or its jet is doubtful. The magnetic 
fields appear to be oriented parallel to the features but due to the range of orientations seen 
in these features, their magnetic fields are not part of a persistent, large scale dipole field. 

Further investigation and observations will be required to understand the details of this 
complex and interesting system. 
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